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Abstract 

We  present  an  asymptotically  derived  boundary  element  method  for  the  Helmholtz  equation  in 
exterior  domains,  in  which  the  basis  functions  are  asymptotically  derived.  Each  basis  function  is  the 
product  of  a  smooth  amplitude  and  an  oscillatory  phase  factor,  like  the  asymptotic  solution.  The  phase 
factor  is  determined  a-priori  by  using  arguments  from  geometrical  optics  and  the  geometrical  theory  of 
diffraction,  while  the  smooth  amplitude  is  represented  by  high  order  splines.  Our  approach  accounts 
for  all  the  components  of  the  scattered  field  namely  the  reflected,  shadow  forming  and  diffracted  fields, 
and  we  demonstrate  that  it  is  substantially  more  accurate  than  an  approach  which  accounts  for  the 
reflected  field  only.  Two  types  of  diffracted  basis  functions  are  presented:  the  first  accounts  for  the 
dominant  oscillatory  behavior  in  the  shadow  region  while  the  second  also  accounts  for  the  decay  of  the 
amplitude  there.  Although  the  method  is  applicable  to  a  variety  of  scatterers,  we  focus  our  attention 
here  on  scattering  from  smooth  convex  bodies  in  two  dimensions.  Our  computations  with  a  conducting 
circular  cylinder  demonstrate  that  the  number  of  unknowns  necessary  to  achieve  a  given  accuracy  with 
this  new  basis  is  virtually  independent  of  the  wave-number. 

1  Introduction 

The  numerical  solution  of  scattering  problems  modeled  by  the  Helmholtz  equation  in  exterior  domains, 
is  exceedingly  difficult  if  not  computationally  intractable  by  standard  computational  schemes,  for  fre¬ 
quencies  in  the  mid  and  high  regime.  Indeed,  for  three  dimensional  problems  the  number  of  unknowns  in 
both  the  finite  difference  and  the  finite  element  methods  scales  like  the  cube  of  ka ,  where  k  is  the  wave 
number  and  a  is  a  typical  dimension  of  the  scatterer,  while  in  a  boundary  integral  method  it  scales  like 
the  square  of  ka.  In  recent  years  various  researchers  [1-10]  have  developed  special  variants  of  the  finite 
element  and  the  boundary  integral  equation  methods  which  reduce  the  computational  complexity  of  the 
problem.  The  common  feature  among  these  algorithms  is  the  use  of  a-priori  knowledge  regarding  the 
oscillatory  nature  of  the  solution  in  the  design  of  the  numerical  scheme. 

In  this  paper  we  develop  a  boundary  element  collocation  method  for  the  integral  formulation  of  the 
Helmholtz  equation,  in  which  the  basis  functions  are  asymptotically  derived.  Each  basis  function  is 
the  product  of  a  smooth  amplitude  and  a  highly  oscillatory  phase  factor,  like  the  asymptotic  solution. 
The  phase  factor  is  determined  a-priori  by  geometrical  optics  (GO)  and  the  the  geometrical  theory 
of  diffraction  (GTD)  [11],  while  the  smooth  amplitude  is  a  high  order  spline.  Each  basis  function  is 
associated  with  a  field:  reflected,  shadow  forming,  diffracted,  multiply  diffracted  etc...  and  we  use 
GO  and  GTD  to  determine  which  fields  are  present.  The  method  is  applicable  to  a  variety  of  scattering 
problems,  but  in  this  paper  we  focus  on  smooth  convex  scatterers  in  two  dimensions,  and  we  demonstrate 
the  method  on  a  conducting  circular  cylinder. 

Our  method  is  analogous  to  the  finite  element  method  which  we  developed  for  the  Helmholtz  equation 
in  [8],  in  which  GO  and  GTD  where  used  to  determine  the  components  of  the  scattered  field,  and  the 
oscillatory  phase  factor  of  the  associated  basis  functions.  It  extends  the  approach  taken  in  [1,  4,  10] 
for  the  integral  equation,  in  that  our  basis  functions  account  for  all  the  fields  present,  including  the 
diffracted  fields.  This  greatly  improves  the  accuracy  of  the  computed  results,  and  enables  the  method 
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to  be  extended  to  problems  where  the  previous  approaches  are  not  effective.  Our  method  is  also  similar 
to  those  presented  in  [2,3,9]  for  discretizing  the  boundary  integral  equation,  with  basis  functions  which 
are  the  product  of  polynomials  and  plane  waves.  The  latter  approach  [9]  draws  it’s  roots  in  the  partition 
of  unity  finite  element  method  [5,12]  for  the  partial  differential  equation,  and  it  makes  use  of  a  multiple 
number  of  plane  waves  on  each  element,  in  directions  which  uniformly  cover  the  unit  circle.  In  our 
approach,  both  the  number  of  basis  functions  and  their  associated  phase  factors  are  determined  by  GO 
and  GTD.  Hence,  the  phase  factors  are  not  necessarily  plane  waves,  and  each  basis  function  is  very 
effective  at  representing  at  least  one  component  of  the  scattered  wave. 


2  Problem  formulation  and  discretization 
2.1  Problem  formulation 

We  consider  scattering  by  a  convex  curve  T.  of  a  time  harmonic  plane  wave  of  unit  amplitude,  exp(ik-x). 
where  k  is  a  vector  of  magnitude  k,  and  x  =  ( x ,  y).  Hence  we  seek  a  solution  to 

A  u  +  k2u  =  0,  (1) 

in  the  unbounded  domain  exterior  to  T.  We  seek  u  of  the  form 

u  =  e*kx  +  tts,  (2) 

where  exp(ik  •  x)  and  us  are  the  incident  and  scattered  fields,  respectively.  Here  we  impose  on  u  a 
Dirichlet  condition  and  on  us  the  Sommerfeld  radiation  condition 


u{x,y)  =  0,  (x,y)  €  T, 


(3) 


and 

limrW  ^-ikA  =  0. 

r-¥  oo  V  or  ) 

The  computational  problem  which  we  wish  to  solve  is  (1)  for  tts,  subject  to  the  conditions 

us  =  —  exp(ik-x),  (x,  y)  E  T, 


(4) 

(5) 


and  (4).  We  derive  an  integral  equation  for  this  problem  in  a  standard  fashion  (see:  [13]), 

/ 


dus(x)  f  <9G(x,£)  exp(ik  •  £) 

v  -G(x,£)dx  =  —  /  exp(dc  ■  x) - ^-^-dx+  s' 

r 


dn 


dn 


(6) 


where  the  free  space  Green’s  function  G(x,  £)  =  (i / 4) Hq  (k\\x  —  Xj||),  and  Hq  is  the  Hankel  function  of 
the  first  kind.  For  the  purpose  of  demonstrating  the  accuracy  of  this  approach,  it  was  not  necessary 
to  consider  more  stable  formulations  such  as  the  combined  fields  approach  [14].  Indeed,  for  the  wave 
numbers  considered  here  the  accuracy  is  excellent  as  we  shall  see. 


2.2  Discretization  by  a  boundary  element  method 

We  shall  approximate  dus  jdn  with  a  boundary  element  method.  Hence,  we  introduce  basis  functions 
Nj(x).  j  =  0, . . . ,  N  —  1  defined  on  T  and  look  for  an  approximation  of  the  form 


dus(x) 

dn 


N- 1 

cjNjix). 

3=o 


(7) 
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In  order  to  obtain  a  system  of  equations  for  the  coefficients  Cj  in  (7)  we  substitute  the  right  side  of  (7)  for 
dus /dn  into  (6),  and  evaluate  the  terms  in  this  equation  at  N  distinct  collocation  points  £o,. . .  ,£n- i- 
This  yields  the  system  of  equations 

Ai,j  =  j  N. j(x)G(x,&)<ix,  9i  =  -  j  exp(ik  •  x)~^~^x  +  6Xp^  ^  .  (8) 

T#  r 

where  T?  is  the  support  of  Nj  on  T. 

We  assume  here  that  the  scatterer  T  is  parametrized  by  the  mapping 

T(s),  a  <  s  <  6,  T(a)  =  T(6).  (9) 

Then,  we  define  the  standard  basis  functions  and  collocation  points  by  first  introducing  the  mesh  points 
Si  i  =  0, . . . ,  N  —  1,  in  [a,  6]  and  the  associated  B-splines  M-(s)  of  polynomial  degree  /,  centered  at  st . 
The  standard  basis  functions  N-.  and  collocation  points  G  are  defined  on  T  by 


&  =  r(ai),  (io) 

JV?(x)  =M/(r"1(x)),  where  x  e  T.  (11) 

2.3  Asymptotically  derived  basis  functions 

The  asymptotic  theory  for  problem  (l)-(4),  based  on  GO  and  the  GTD,  is  described  in  detail  in  [15,16]. 
Here,  we  review  the  essential  elements  of  that  theory  necessary  for  the  development  of  asymptotic  basis 
functions.  Additional  relevant  details  can  be  found  in  [8, 17].  In  the  asymptotic  theory,  one  seeks  the 
solution  us  in  (2)  as  a  superposition  of  the  form 

L 

us{x,y)  =  Ai(x,y)exp(ikSi(x,y)).  (12) 

i=i 

Here  L  is  the  number  of  fields,  and  for  a  smooth  convex  scatterer  there  are  four  terms 

us  =  ur  +  usf  +  ud+  +  ud~ .  (13) 

They  are  the  reflected  field  tir,  the  shadow  forming  field  ,  and  two  diffracted  fields  ud+ ,  ud~  which 
are  associated  with  special  points  of  T.  Specifically,  each  incident  ray  which  is  tangent  to  the  cylinder 
gives  rise  to  a  surface  diffracted  ray  which  travels  around  T  starting  at  the  point  of  diffraction.  When  the 
scatterer  is  a  circular  cylinder  of  radius  1  centered  at  (0,0)  and  k  =  (1,0)  in  (2),  the  points  of  tangency 
are  (0, 1)  and  (0,  —1),  respectively.  Moreover,  usf  satisfies  usf  +  ul  =  0  in  the  region  x  >  0,  \y\  <  1,  so 
usf  =  —  exp (ikx)  there.  In  our  formulations  it  contributes  an  additional  term 

—k/ 4  J  cos #exp(?£;x)i7o(fc||x  —  x;||)dx,  9  =  arg(x), 

r\rr 


to  gi  in  (8). 

Each  field  in  (13)  has  the  form 


A(x,y)exp{ikS{x1y))1 


(14) 
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where  the  amplitude  A  is  smooth  and  the  phase  S(x,y)  can  be  determined  by  the  ray  method.  Only  the 
reflected  and  shadow  forming  fields  are  non-zero  on  the  surface  of  the  scatterer  but  the  normal  derivative 
of  all  fields  is  non-zero  there.  With  T  the  unit  circle 

T(s)  =  (coss,sins),  (15) 

and  k  =  (1,0)  in  (2)  we  have 

Sr{x,y)  =  x,  Ssf{x,y)  =  x,  Sd+  =  s+{x,y),  Sd  =s~(x,y),  s+,s~>  0.  (16) 

The  functions  s+  and  s~  are  respectively  the  distances  traveled  around  the  cylinder  by  the  surface 
diffracted  rays  from  the  diffraction  points  (0, 1)  and  (0,  —1)  in  the  clockwise  and  counterclockwise  direc¬ 
tions  to  the  point  (x,  y).  The  functions  .s+,  s~  are  multivalued  and  unbounded  because  surface  diffracted 
rays  travel  an  infinite  number  of  times  around  the  cylinder.  This  is  easily  overcome  as  explained  in  [8, 17], 
and  .s+,  s~  in  (16)  are  the  shortest  distances.  Moreover,  the  quantities  dud+  /dn  and  dud  /dn  decay 
exponentially  by  a  damping  factor  D(s+)  and  D(s~),  respectively  (see:  [16,17])  where 

D(s )  =  exp(iroA:1'/3s),  and  to  =  exp(*7r/3.0)  (9tt/2)2^3  /4.O.  (17) 

It  follows  that  dud+  /dn  is  exponentially  smaller  than  dud  /dn  when  s+  >  7t/2  and  k  — »  oo,  and  vice 
versa  when  s~  >  7t/2.  However,  both  dud+  /dn  and  dud  /dn  will  be  non-negligible  compared  to  each 
other  in  the  neighborhood  of  (x,y)  =  (1,0)  where  s+  =  s~  =  tt/2. 

This  leads  to  a  decomposition  of  T  into  four  zones 

rr  =  (T(s)  |  7t/2  <  s  <  37t/2}  ,  rd+  =  (r(s)  I  a  <  s  <  vr/2}  ,  (18) 

rd  =  (r(a)  |  37r/2  <  s  <  2ir  —  a}  ,  T0  =  (T(s)  |  —a  <  s  <  a}  ,  (19) 

with  a  a  small  positive  scalar.  In  each  zone  we  use  basis  functions  of  the  form 

Ni(j,  l,  t,  x,  k)  =  IVj(x)  exp(?A:5't(x)),  for  x  £  T(  te{d+,d_,r},  (20) 

where  5*  is  defined  in  (16),  A^(x)  in  (11)  and  in  (18)  and  (19).  In  the  overlap  region  T0  we  use 
basis  functions  associated  with  both  Td+  and  Td-.  We  also  introduce  the  modified  basis  functions  in 
Td+,  T d-,  T0,  which  include  the  decay  (17) 

^(j, /,  t,  x,  k)  =  ^Vj(x)  exp(ikSt(x.))D(St(x.)),  for  x  E  t€{d+,d~).  (21) 

In  our  calculations  we  shall  compare  the  performance  of  the  basis  functions  (20)  and  (21)  with 

N3(j,l,x,k)=Nj(x)exp(ikx),  x  e  T,  (22) 

which  are  based  on  the  GO  ansatz  used  in  [1,4,10]. 

2.4  Implementation  details 

Several  issues  need  to  be  considered  in  a  numerical  implementation  in  order  for  the  method  to  be  effi¬ 
cient  and  accurate,  and  they  are  presented  in  detail  in  [17].  The  numerical  integration  of  the  oscillatory 
integrals  in  (8)  are  computed  by  specialized  quadrature  rules,  in  time  virtually  independent  of  k.  Precon¬ 
ditioning  of  the  system  (8)  by  the  norm  of  each  column,  when  basis  functions  (21)  are  used,  is  required. 
Improved  efficiency  can  be  achieved  by  mesh  adaptivity,  as  discussed  in  [8],  and  the  use  of  variable  or¬ 
der  splines.  Appropriate  measures  need  to  be  taken  to  avoid  dependence  between  basis  function  in  the 
overlap  region  for  small  values  of  &,  as  described  in  [8, 17]. 
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3  Numerical  calculations 

d(us — US  f } 

We  now  compare  the  accuracy  of  the  numerical  approximation  to  v  dn — -  obtained  with  each  of  the 
three  families  of  basis  functions  (20),  (21)  and  (22),  in  which  JVj  are  B-splines  of  order  l  =  5.  The  number 
of  discretization  points  is  250. 

Figure  1  presents  the  absolute  value  of  the  real  and  imaginary  parts  of  the  exact  and  computed 
solutions  at  the  nodes  of  the  grid  in  the  shadow  region,  for  k  =  750.  In  most  of  the  shadow  region  bases 
(20)  and  (21)  yield  excellent  results  while  basis  (22)  provides  an  accurate  solution  in  only  half  of  that 
region.  Basis  (21)  is  slightly  better  than  (20).  All  three  approximations  agree  to  graphical  accuracy  in 
the  illuminated  region  which  we  do  not  show  here.  Basis  (22)  is  accurate  near  the  shadow  boundary 
because  the  phase  of  the  diffrcated  field,  s+,  is  well  approximated  there  by  x ,  as  indicated  by  the  relation 

+3 

x  =  cos(-7r/2  —  s+)  =  s+  —  ^ — I-  0(s+5).  (23) 


Angle/n  Angle/n 


Figure  1:  The  absolute  value  of  the  real  and  imaginary  parts  of  the  exact  and  computed  values  of 
d  (us  —  usf)  /dn  for  k  =  750,  in  the  left  and  right  panels  respectively.  The  exact  solution,  and  approxi¬ 
mations  with  basis  (20)-(22)  are  plotted  with  o,  +,  *,  and  x,  respectively. 
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Table  2  presents  the  error  between  exact  and  computed  values  of  v  dn — -,  on  the  cylinder  surface 
over  the  range  k  =  10  —  5000,  evaluated  at  the  nodes  of  the  grid.  Both  the  maximum  and  median  errors 
are  presented.  The  superiority  of  bases  (20)  and  (21)  over  (22)  is  apparent.  Basis  (21)  provides  an  order 
of  magnitude  greater  accuracy  than  basis  (21). 


k 

Basis  (20) 

error  oo 
Basis  (21) 

Basis  (22) 

i 

Basis  (20) 

/ledian  |  error 
Basis  (21) 

Basis  (22) 

10 

2.66e-07 

4.82e-07 

1.19e-08 

3.45e-09 

6.88e-09 

5.52e-10 

50 

7.67e-04 

9.88e-04 

8.54e-06 

6.02e-05 

7.71e-05 

5.30e-07 

100 

6.53e-04 

5.61e-04 

3.07e-04 

3.47e-05 

3.44e-05 

3.20e-05 

250 

6.99e-05 

3.82e-05 

3.42e-02 

3.91e-06 

1.60e-06 

4.26e-03 

500 

1.68e-05 

1.73e-06 

5.63e-02 

6.96e-07 

9.57e-08 

5.83e-03 

1000 

1.09e-04 

1.49e-05 

5.41e-01 

3.06e-06 

6.27e-07 

1.01e-01 

2000 

7.20e-04 

1.34e-04 

3.99e+00 

1.49e-05 

4.37e-06 

9.29e-01 

2500 

1.33e-03 

2.63e-04 

1.16e+01 

2.33e-05 

9.81e-06 

2.12e+00 

5000 

8.99e-03 

2.24e-03 

2.31e+01 

4.30e-04 

7.13e-05 

4.83e+00 

Figure  2:  The  maximum  and  median  error  at  the  discretization  points  on  the  surface  of  the  cylinder. 

The  increase  in  the  error  as  a  function  of  k  which  can  be  seen  in  Table  2,  is  due  to  the  fact  that 

,  is  proportional  to  k  in  the  illuminated  region.  This  does  not  affect  by  much  the  error  in 

the  total  field  u1  +  us  in  the  exterior  of  the  cylinder  as  indicated  in  Table  3.  This  table  presents  the 
maximum  and  median  relative  error  in  | u1  +  us\  for  a  set  of  points  A  in  the  illuminated  region  exterior 
to  the  cylinder  defined  by 

A  =  { (ri  cos  9j ,  Ti  sin  9j)  \  n  =  1  +  .5  i,i  =  1,  -  ,4,  6j  =  j2n/17}\{(x,y)  \  x  <  0,  \y\  <  1}.  (24) 

Here  again  we  see  that  bases  (20)  and  (21)  are  at  least  4  orders  of  magnitude  more  accurate  than  the 
basis  (22).  Tables  3  presents  the  maximum  and  median  relative  error  in  \ul  +  us\  for  a  set  of  points  B  in 
the  shaded  region  exterior  to  the  cylinder  defined  by 

B  =  {(ri  cos  ^  sin 9j)  \  rt  =  1  +  .5?,  i  =  1,  •  •  •  ,4,  9j  =  j27r/17 }  fl  {(rr,  y)  \  x  <  3,  \y\  <  1}.  (25) 

We  excluded  from  B  points  in  the  deep  shadow  close  to  the  cylinder  surface,  as  we  did  not  expect  basis 
(22)  to  perform  well  there.  However,  we  see  that  throughout  the  shadow  region  it  performs  quite  poorly, 
Basis  (21)  on  the  other  hand  is  accurate  to  at  least  two  significant  figures  even  for  k  =  5000.  The 
reduction  in  the  number  of  significant  digits  in  the  shadow  region  as  k  increases  is  explained  in  [17]. 


d(us-usf) 

dn 
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k 

relative  error 

OO 

Median  |relative  error| 

Basis  (20) 

Basis  (21) 

Basis  (22) 

Basis  (20) 

Basis  (21) 

Basis  (22) 

1.96e-10 

4.02e-10 

1.41e-ll 

7.99e-12 

1.63e-ll 

1.18e-12 

1.80e-06 

2.48e-06 

1.36e-08 

6.73e-08 

l.lle-07 

1.40e-09 

1.14e-06 

l.lle-06 

2.70e-06 

1.22e-07 

1.23e-07 

4.16e-07 

5.26e-08 

2.94e-08 

4.76e-05 

6.26e-09 

3.25e-09 

5.98e-06 

9.33e-09 

1.18e-09 

8.37e-05 

9.08e-10 

7.52e-ll 

5.21e-06 

1000 

2.06e-08 

3.96e-09 

5.49e-04 

2.43e-09 

5.24e-10 

4.99e-05 

2000 

9.56e-08 

1.31e-08 

1.44e-03 

2.27e-09 

4.82e-10 

1.87e-04 

2500 

1.91e-07 

2.45e-08 

1.95e-03 

8.69e-09 

1.89e-09 

3.82e-04 

5000 

7.17e-07 

8.29e-08 

2.59e-03 

1.87e-08 

3.29e-09 

4.33e-04 

Figure  3:  The  maximum  and  median  relative  error  on  the  set  A  of  points  in  the  illuminated  region. 


k 

relative  error 

OO 

Median  |relative  error| 

Basis  (20) 

Basis  (21) 

Basis  (22) 

Basis  (20) 

Basis  (21) 

Basis  (22) 

4.51e-10 

9.13e-10 

1.14e-ll 

2.13e-10 

4.31e-10 

4.46e-12 

3.99e-06 

4.56e-06 

2.37e-08 

1.40e-06 

1.60e-06 

9.68e-09 

1.30e-05 

1.68e-05 

1.70e-05 

2.24e-06 

2.69e-06 

3.00e-06 

5.65e-06 

3.05e-06 

2.99e-03 

8.03e-07 

3.29e-07 

3.93e-04 

2.15e-07 

3.05e-08 

2.89e-02 

4.99e-08 

9.51e-09 

2.77e-03 

1000 

8.21e-06 

2.43e-06 

4.55e-01 

7.33e-07 

1.77e-07 

4.28e-02 

2000 

7.86e-04 

1.37e-04 

8.25e+00 

4.50e-06 

1.14e-06 

2.57e-01 

2500 

4.01e-03 

4.73e-04 

7.98e+01 

3.72e-05 

9.03e-06 

1.03e+00 

5000 

1.80e-01 

1.61e-02 

2.11e+02 

4.44e-04 

1.25e-04 

3.28e+00 

Figure  4:  The  maximum  and  median  relative  error  on  the  set  B  of  points  in  the  shaded  region. 
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